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This technical report presents an algorithm and code to solve for two similar 
classes of problems. Given a list of points, find the equations for the best line 
that passes a minimum distance from each point. Given a list of lines, find 
the point that is a minimum distance from each line. Singular value decom- 
position is used to solve for the inconsistant matrix. Numerical robustness is 
tested for alternate solutions to the line problem, and the best solution form 
is returned. This avoids singularities in best fit lines that are oriented along 
certain directions. This algorithm is implemented in “C”, and has applications 
in locating points in tomographic slice stacks, and finding trajectories of high 
energy particles passing through arrays of detectors. 

The source code is presented in the appendix. 


• Find the line that has the minimum RMS 1 distance from a list of points 
denoted (a?,- , y,- , z,) , » = 1, 2, 3, . . . 

• Find that point that has the mimimum RMS distance from a given list of 
lines. 

*This work was supported by an equipment loan from Silicon Graphics, Inc. D. Karron is 
with Department of Applied Science at New York University and the Department of Surgery 
at New York University Medical Center, 560 First Avenue, New York, New York, 10016. 
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Abstract 


An algorithm is presented to solve for two similar classes of problems. 
Given a list of points, find the equations for the best line that passes a 
minimum distance from each point. Given a list of lines, find the point 
that is a minimum distance from each line. Singular value decomposition 
is used to solve for the inconsistant matrix. Numerical robustness is tested 
for alternate solutions to the line problem, and the best solution form 
is returned. This avoids singularities in best fit lines that are oriented 
along certain directions. This algorithm is implemented in “C” , and has 
applications in locating points in tomographic slice stacks, and finding 
trajectories of high energy particles passing through arrays of detectors. 
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1 Source Code for Subroutines 


/* This code a nd the underlying algorithm were written and developed by 

* Dan Karron (karron@nyu.edu). You are free to use it, but please 

* acknowledge me if you do anything worth while with it. If it dosent 

* work for you, let me know and I will try to help. It works quite well for 

* me and my colleagues here. July 1991. 

*/ 


#include "std r,i 
^include "stdl r,; 

^include "s r-i 
#include "gl.h" 

^include "str rH 
^include "math.h" 

#include "s r,i 

#include "nr .h" /* requires the Numerical Recipies library in C */ 

^include "nmt^ 

^include "BestL r,i 
#include "Pr- h" 

#define SVD_EDIT 


)**********************************************************************^************/ 
static double RunCalculation(float **a,int m, float *x) 

{ 

/* Calculate a best solution, then return with a figure of merit for 

* how good the solution really is by using the solution and taking the 

* difference between the solution and the known arbitrary line point normal 

* constant 1.0 . 

*/ 


register int ij; 

static int n=2; /* solve in 2d subspace */ 
double wmax,wmin, condition; 
float **u,*w,**v,*b,*c; 


u=matrix( 1 ,m, l,n) ; 
v=matrix( 1 ,m, 1 ,n) ; 
w= vector (l,m); 
b=vector(l,m); 
c=vector(l,m); 


for(i=l;i<m;i++) 
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It is crucial to represent each line as the intersection of two planes. The 
solution of the second problem then follows simply from the results of the first 
problem. Each line is given as a pair of intersecting plane equations : 

A Xi + B yi + C Z{ = 1 

D X{ + E yi + F Zi = 1 

where » = 1, 2, 3, . . . , n 

The first problem can be formulated as 


[A] [X] = [B] 


where 
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Note that the data matrix [j 4] has, in general, many more rows than columns, 
so that this system is usually over determined. The SVD (Singular Value De- 
composition) solution is thus one of many possible solutions and can be shown 
to be the one with the required RMS error properties. 

The SVD technique is based on the observation that any matrix [.A] whose 
number of rows N is greater or equal to its number of columns N , can be 
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diagonal matrix [W], with positive or 


written as the product of an N x M column orthogonal matrix [Z7], an M x M 


zero elements, and the transpose of an 


M x M orthogonal matrix [V] [1], [2], M3]; that is, 


[A] = mw][v\ T . 


The “standard” SVD solution, which can be verified by back substitution, 
is given as : 


m = mw W[s] 


The problem with this solution is that the pairs of planes that result can 
be close to parallel, leading to severe numerical instabilities if these values are 
used to solve the second problem. To minimize these numerical instabilities we 
constrain the solution planes to be perpendicular. Any calculation with these 
lines will suffer from numerical instability if the rows of the line matrix are not 
“maximally” linearly independent. 

Intuatively, think of the line as two intersecting planes. Think of each plane 
as minimally defined as three points. Now at first glance it would appear that 
6 points are required to define two planes. Remember that the planes can have 
points in common, and that we can take two triangles for two planes that share 
two points. So we can minimally solve for two planes in space with 4 points. 
If we now constrain our planes to be perpendicular, then we do not need one 
point. Now if we also allow our perpendicular planes to rotate about the line, 
then we do not need the remaining third point to constrain the perpendicular 
planes in any orientation with respect to our line. Therefore, we can use two 
points to define a line, but the orientations of the planes are undefined, but 
perpendicular to each other. 

From the above argument, it is clear that there are multiple equivalent pairs 
of planar equations for any line. Perpendicularity means that any one coefficient 
from the first plane {A, J5, C} and any one coefficient not in the same column 
{ D , E , F} can be 0. This is equivalent to projecting the line onto any two of the 
x — y,y — ZyOtz — x planes. This reduces the 2xn overdetermined equations in 
6 unknowns to a three systems of n overdetermined equations in two unknowns. 
However, if any of the z,- ,y,- , or z, values for our lines approach zero, then 
one pair of equations becomes singular. To avoid singularity in any particular 
direction, we must determine which three pairs of equations to use. Since each 
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equation solution must be equal to 1, any difference from 1 is numerical insta- 
bility symptomatic of singularity. We throw out the equation pair with largest 
error, and keep the two pairs of equations with smaller error. A line that has 
a direction that approaches singularity will have an bigger numerical error in a 
trial solution than the others. However, we do not know a priori which of the 
systems of equations to use. We try all three systems and throw out the system 
with the largest error. 


The solution of the second problem is given by 

[A]X = B 


Where 




Aq Bo Co 
Dq Eq Fo 
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[B] 


X = 


This identifies a single point that best solves the lines by minimizing the 
distance to each of the lines. The point can be interperted as the singular value 
decomposition solution of the best solution to the system of planes describing 
the lines. 
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* coefficients (where the point constant is taken as 1.0, always!) 

* in lcoeff for two point normal planes that describe the line that 

* best fits the list of points in three dimensions. 

V 

register int i; 
float **a,*x[4]; 
double c[4]; 
int cmaxi; 
double cmaxf; 
static int n=2; 

#ifdef DEBUG 
Debug=2; 

#endif DEBUG 

#ifdef DEBUG 

for(i=0;i<m;i++) /* Recite the incomming point list just to be certain. *f 

{ 

printf("NueITogg r ' i nocker Po N ,*/,! ,y,f\n",i> 

*(plist[i]), *(plist[i]+l), *(plist[i]+2)); 

} 

#endif DEBUG 

if( m < n ) 

{ 

fprintf(stderr, H too lew po r<i '/,d<'/,d to do anytlr* M<N\n",m,n); 
return; 

} 

a=matrix( 1 ,m, 1 ,n) ; 

x[l]=vector(l,n); /* hang three solution vectors off x */ 

x[2]= vector (l,n); 

x[3]=vector(l,n); 

/♦ first try */ 
for(i=0;i<m;i++) 

{ 

a[i+l][!]= *(plist[ij); /♦ x */ 

a[i+l][2]= *(plist[i]+l);/f y */ 

} 


#ifdef DEBUG 
printf("l r-i try\n M ); 


rs 
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{ . 

b[i]=1.0; /* point normal plane constant */ 
for(j=l j<2y++) 

{ 

u [i] [i] — a (i] [j] ; A don’t touch matrix a */ 

} 

} 

#ifdef DEBUG 
/* Lets see the intput */ 

Print_Matrix(" r,i M ,l,a,l,m,l,n); 

Print_Vector(" constant s " , l,b, 1 ,m) ; 

#endif DEBUG 

svdcmp(u l m,n l w > v); /* decompose u into u,v, and w */ 

wmax=0.0; /* edit diagnonal which is vector v */ 
for (i= 1 ;i<m;i++ ) 

{ 

if(w[i] > wmax) 

{ 

wmax=w[i]; 

} 

} 

#ifdef DEBUG 

printf ( " wmax='/,l \n" , wmax) ; 

#endif DEBUG 

wmin= wmax* ( 1 . Oe — 6 ) ; 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 

{ 

w[i]=0.0; 

} 

} 

#ifdef DEBUG 
printf("wm r,i n",wmin); 

#endif DEBUG 


svbksb^^jVjm^jbjX); f* svd back substitute to solve for best x */ 
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#ifdef DEBUG 

Print_Matrix( " lbs " , 1 ,u, 1 , m, 1 ,n) ; 

Print_Vector("d r - i ,l,w,l,m); 

Print_Matrix( "rhs " , 1 ,v, l,m, 1 ,n) ; 

Print_Vector("solut r,i ,l,x,l,n); 

#endif DEBUG 

nrdot (a, 1 ,m f 1 ,n,x,c) ; 

#ifdef DEBUG 

Print_Vector("test solut^ )",l,c f l,m); 

#endif DEBUG 

condition=0.0; /* this is not really condition, just an figure of merit 
* for the absolute error, as absolute difference from 1.0 

V 

for(i=l;i<m;i++) 

{ 

if((c[i]— 1.0)<0.0) 

condition — = (c[i] — 1.0); 
else 

condition += (c[i] — 1.0); 

} 

#ifdef DEBUG 

printf( " cond r ^ n=*/,l \n M , condition) ; 

#endif DEBUG 


free_matrix(u, 1 ,m, 1 ,n) ; 
free_matrix(v, 1 ,m, 1 ,n) ; 
free_vector(b, 1 ,m) ; 
free-vector(c,l,m); 
free_vector( w, 1 ,n) ; 


return condition; 

} 

y^*********:M**********************************************************************y 

void NueNogginKnocker(int m,float **plist,float lcoeff[2][3]) 

{ 

A 

* Given a list of m point triples (x,y,z) ( which actually should be declared 

* (*P)ft] blit I had trouble making that work) return the point normal line 
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free_vector(x[2] , 1 ,n) ; 
free_vector(x[3] , 1 ,n) ; 

} 

^*********************************************************************************y 


#ifdef OBSOLETE 

/* this version does not check for best two solutions out of three, 
* and has a bad numerical instability in certain directions. 

*/ 


void NogginKnocker(int m, float **plist, float lcoeff[2][3]) 
{/* project onto two planes, first pass= x,y *j 
/* second pass, x,z *j 
register int ij,k; 

float **u l *w,**v l **a,*b,*c,*x,**cvm; 
float wmax,wmin; 
int n=2; 

float al ,a2,bl,c2; 


#ifdef DEBUG 
Debug=0; 

#endif DEBUG 

for(j=0;j<m;j++) 

{ 

printf("Po r ' */,f , */,f , 

*(plist[j]), *(plist£j]+l), *(plist[j]+2)); 

} 

if(m < n ) 

{ 

fprintf(stderr, M Umlout M<N\n H ); 
return; 

} 

a=matrix( 1 ,m, 1 ,n) ; 
b= vector (l,m); 
c=vector(l,m); 
u=matrix( 1 ,m, 1 ,n) ; 
v=matrix(l,m,l,n); 
cvm=matrix( 1 ,m, 1 ,n) ; 
w=vector(l,m); 
x=vector(l,n); 


for(i=0;i<m;i++) 
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#endif DEBUG 

c[l]=RunCalculation(a,m,x[l]); 

/* second try */ 

for(i=0;i<m;i++) 

{ 

a[i+l][l]= *(plist[i]); /* x */ 
a[i+l][2]= *(plist[i]+2);/* z */ 

} 

#ifdef DEBUG 
printf(" second try\n"); 

#endif DEBUG 

c [2] =RunCalculation(a l m J x[2] ) ; 

/* third try */ 

for(i=0;i<m;i++) 

{ 

a[i+l][l]= *(plist[i]+l); /* y */ 
a[i+l][2]= *(plist[i]+2); /* z */ 

} 

#ifdef DEBUG 

printf("th r ^ try\n"); r 

#endif DEBUG 

c[3]=RunCalculation(a,m,x[3]); 

cmaxf=0.0; /* which is the worst of the three ?, keep the two best ! */ 
for(i=l;i<3;i++) 

{ 

if(cmaxf<c[i]) 

{ 

cmaxf=c[i]; 

cmaxi=i; 

} 

} 

#ifdef DEBUG 

printf("cmax=’/,f , cmax r,i ='/.d\n H ,cinaxf,cmaxi); 

#endif DEBUG 
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switch(cmaxi) /* wire the best two solutions 

* back into the solution for the line 

*/ 

{ 

case 3: 

lcoeff[0] [0] =x[l][l] ; 
lcoeff[0] [1] =x[l] [2] ; 
lcoeff[0][2]=0.0; 

lcoeff[l] [0] =x[2][l] ; 
lcoeff[l][l]=0.0; 
lcoeff[l] [2] =x[2][2] ; 
break; 

case 2: 

lcoeff[0] [0] =x[l][l] ; 
lcoeff[0][l]=x[l][2]; 
lcoeff[0][2]=0.0; 

lcoefF[l][0]=0.0; 
lcoeff[lj[l]=x[3][l]; 
lcoefF[l] [2] =x[3] [2] ; 
break; 

case 1: 

lcoeff[0][0]=x[2][l]; 
lcoeff[0][l]=0.0; 
lcoefF[0] [2] =x [2] [2] ; 

lcoefF[l][0]=0.0; 
lcoeff[l][l]=x[3][l]; 
lcoeff[l] [2] =x[3] [2] ; 
break; 

} 

#ifdef DEBUG 

printf("lcoe« [0] = 7.1 ,7.1 , 7.1\n" ,lcoeff[0] [0] ,lcoeff[0] [l],lcoeff[0] [2] ) ; 
printf("lcoeif [1] = 7.1 ,7.1 ,y.i\n H ,lcoefF[l][0],lcoefF[l][l],lcoeff[l][2j); 
#endif DEBUG 


free_matr ix(a, 1 ,m, 1 ,n) ; 
free_vector(x[l] , 1 ,n) ; 
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free_vector(c,l,m); 
free_matrix(u, 1 ,m, 1 ,n); 
free_matrix(v, 1 ,m, 1 ,n) ; 
free_vector( w, 1 ,n) ; 
free_vector(x, 1 ,n); 
free_matrix(c vm, 1 ,m, 1 , n) ; 

#ifdef DEBUG 
Debug=0; 

#endif DEBUG 


lcoeff[0][0]=al; 

lcoeff[0][l]=bl; 

lcoeff[0][2]=0.0; 

lcoeff[l][0]=a2; 

lcoeff[l][l]=0.0; 

lcoefF[l][2]=c2; 

} 

#endif OBSOLETE 


/**************************:^*************************************************y 

void NitherNoid(int mm,float **llist, float point[3],int debug_me) 

{ 

/* Find the best point to solve a mm long list of lines attached to llist. 

* Actually, the solution is defined the same as the best point to solve 

* for a list of lines consisting of pairs of planes. 

*/ 

register int ij,k; 

float **u,*w l **v J **a l *b,*c > *x,**cvm; 
float wmax,wmin; 

int n=3; /* solution is in 3 dimensions */ 

int m=mm*2; /* each line consists of two planes */ 


#ifdef DEBUG 
Debug=l; 
#endif DEBUG 


#ifdef DEBUG 

if(debug_me) 

forQ=Oy<mmij++) 

{ 

printf("N r < :L r ’ e['/,d] :*/.+e, , /,+e,y,+e\n\t'/,+e,'/.+e, , /,+e\n"J, 

*(llistp]+0),*(llistp]+l) J *(llist[j]+2), 
*(llist[j]-h3),*(llist[3]-f-4),*(llist[j3-J-5)); 


} 
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{ 

b[i+l]=1.0; 

u [i+l][l]= a [i+l][l]= *(plist[i]); /* x */ 
u [i+ 1] [2]=a[i+ 1] [2] = *(plist[i]+l);A y */ 

} 

Print_Matrix("F r,i Pass, r ’ ,l,a,l,m,l,n); 

Print_Vector("F r,i Pass, constants*', l,b,l,m); 

svdcmp(u,m,n J w,v); 

wmax=0.0; 

for(i=l;i<m;i++) 

{ 

if(w[i] > wmax) 

{ 

wmax=w[i]; 

} 

} 

#ifdef DEBUG 

printf("F r,i Pass wmax=y,f \n",wmax); rs 

#endif DEBUG 
wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 

{ 

w[i]=0.0; 

} 

} 

svbksb(u,w f v,m,n J b,x); 

#ifdef DEBUG 

Print_Matrix( " lbs " , 1 ,u, 1 ,m, 1 ,n) ; 

Print_Vector("d rs 
Print_Matrix( " rhs " , 1 , v, 1 ,m, 1 ,n) ; 

Print_Vector("F^ pass solut^ ,l,x,l,n); 
nrdot(a,l,m,l,n,x,c); 

Print-Vector^'F^ pass test solut” ?)",l,c,l,m); 

#endif DEBUG 

al=x[l]; 

bl=x[2j; 


/* second pass */ 
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} 

} 

wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 

{ 

w[i]=0.0; 

} 

} 

#endif SVDJEDIT 


svbksb(u,w,v,m,n,b,x); 


#ifdef DEBUG 

Print_Matrix( H lh.s ,, J l l u J l l m,l,n); 

Print-Vector^'d” jljWjl.m); 

Print_Matrix( "rhs " , 1 ,v, l,m, l,n); 

Print_Vector("N” solut” ,l,x, l,n); 

nrdot(a,l,m,l,n,x,c); 

Print_Vector("N” test solut” ?) ",1,0,1, m); 

#endif DEBUG 


point[0]=x[l]; 

point[l]=x[2]; 

point[2]=x[3]; 


free_matrix(a,l,m,l,n); 
free_vector(b,l,m); 
free_vector(c , 1 ,m) ; 
free_matrix(u,l,m,l,n); 
free_matrix(v,l,m,l,n); 
free_vector( w, 1 ,n) ; 
free_vector(x, 1 ,n) ; 
free_matrix(c vm, 1 ,m, 1 ,n) ; 
} 




References 
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for(i=0;i<m;i++) 

{ 

u[i+l][l]=a[i+l][l]= *(plist[ij); /* x */ 
u[i+l][2]=a[i-f 1][2]= *(plist[i]+2);/* z */ 

} 

Print_Matrix(" Second Pass, r<i ,l,a,l,m,l,n); 

Print_Vector("Second Pass, cons t ant s",l,b,l,m); 
svdcmp(u,m,n,w,v); 
w max =0.0; 
for(i=l;i<m;i++) 

{ 

if(w[i] > wmax) 

{ 

wmax=w[i]; 

} 

} 

#ifdef DEBUG 

printf("Secpmd Pass wmax='/,f \n",wmax); 

#endif DEBUG 
wmin=wmax*(1.0e— 6); 
for(i=l;i<m;i++) 

{ 

if(w[i] < wmin) 

{ 

w[i]=0.0; 

} 

} 

svbksb^jW.VjmjnjbjX); 

#ifdef DEBUG 

Print_Matrix( " lhs " , 1 ,u, 1 ,m, 1 ,n) ; 

Print_Vector("d r,i ,l,w,l,m); 

Print_Matrix( "rhs ", 1 ,v, l,m, 1 ,n) ; 

Print_Vector("Second pass solnt r,i ,l,x,l,n); 
nrdot(a,l,m,l,n,x,c); 

Print_Vector("Second pass test solut^ b?)",l,c,l,m); 

#endif DEBUG 


a2=x[l]; 

c2=x[2]; 


free_matrix(a, 1 ,m, 1 ,n) ; 
free_vector(b,l,m); 


Karron 


Page 16 


#endif DEBUG 

if(m < n ) 

{ # 

fprintf(stderr,"Too lew plane pa r<i to solve lor M<N\n")> 
return; 

} 

a=matrix(l,m,l,n); 
b=vector(l,m); 
c= vector (l,m); 
u=matrix( 1 ,m, l,n); 
v=matrix(l J m > l,n); 
c vm=matrix( 1 jin, 1 ,n) ; 
w=vector(l,m); 
x=vector(l,n); 

for(i=0;i<mm;i++) 

{ 

b[i*2+l]=b[i*2+2]=1.0; /* point normal plane constant */ 
u[i*2+lj[l]=a[i*2+l][lj= *(llist[i]+0); 

u[i*2-fl][2]=a[i*2+l][2]= *(llist[i]+l); 
u[i*2+l][3]=a[i*2+l][3]= *(llist[i]+2); 

u[i*2+2][l]=a[i*2+2][lj= *(llist[i]+3); 
u[i*2+2][2]=a[i*2+2][2]= * (Hist [i] 4-4); 
u[i*24-2][3]=a[i*24-2][3]= *(llist[i]4-5); 

} 

#ifdef DEBUG 

Print_Matrix("N r,i r<i put",l,a,l,m,l,n); 

Print_Vector( M N r,i s t ant s ” ,1 ,b, 1 ,m); 

#endif DEBUG 

svdcmp^^jnjWjv); 

#ifdefSVD_EDIT 
w max =0.0; 
for(i=l;i<m;i4-4-) 

{ 

if(w[i] > wmax) 

{ 

wmax=w[i]; 
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1 Perspective 

As the title of this article suggests, studies of the magnetic field of the human brain have 
provided a new and important method for “seeing” neuronal activity. These studies make it 
possible to locate specific brain functions and define some of their key functional properties 
as well. However, the title is a double-entendre, for we also have chosen to illustrate the 
power of this technique with examples taken from our group’s studies of the human visual 
system. Worldwide, more attention has been devoted to studying vision than to any other 
sensory modality. On the other hand, researchers in neuromagnetism have devoted more 
attention to the auditory system. Perhaps this is because it is accessible for study with 
stimuli provided in a comparatively simple way (with earphone or speaker). By comparison, 
many forms of visual stimulation (e.g., a computer screen) may impose magnetic artifacts 
that are difficult to avoid. This imbalance in research effort is beginning to change, and visual 
studies have become one of the prime interests of Professor Olli Lounasmaa. His group in 
neuromagnetism recently began research in this area at the Low Temperature Laboratory. 
For this reason, our choice of topic is particularly appropriate for this volume. 

We shall focus on properties of the visual cortex, where signals are sent from the retina 
after first being processed in the thalamus. Only a small part of visual cortex extends across 
the outer surface of the brain, near the most posterior point of the head. The greater part 
continues inward along the left and right walls of the longitudinal fissure, which divides the 
two hemispheres, and reaches about half-way to the center of the head. Looking at the brain 
from the back, the visual areas of left and right banks of the longitudinal fissure fold into the 
respective hemispheres, forming left and right calcarine fissures. The vertical longitudinal 
fissure and horizontal left and right calcarine fissures form a cross-like pattern, so that the 
geometry of visual cortex is often described by a “cruciform model”. 

2 Earliest Measurements 

The first detection of an evoked field arising from the bcain was achieved by Brenner et al. 
in 1974 using visual stimuli. [1] We used a two-hole rf SQUID of Zimmerman design and a 
home-made second-order gradiometer as the detection coil. These measurements were done 
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in a normal laboratory at NYU without a magnetic shield. Additional balancing coils were 
connected in series with the detection coil to improve rejection of environmental field noise. 
The object of our study was to investigate the spatial pattern of magnetic field that appears 
across the scalp when a grating pattern is shown to a subject. 

Because of the high instrumentation noise level at that time (typically 75/T Hz 1 / 2 , one 
order of magnitude greater than today!), the steady-state method was used to achieve a very 
narrow measurement bandwidth. In the steady-state method a given stimulus is presented 
repetitively and the waveform of the response is averaged across individual epochs. If stimuli 
are presented rapidly, say at 5 Hz or greater, the response to one stimulus has not been 
completed before the next stimulus is presented. The strongest Fourier component of the 
averaged waveform has this fundamental period, which can be measured with a lock-in 
am plifier (phase sensitive detector) or multichannel averager. With a lock-in, the bandwidth 
of the meas ure men t is determined by the time constant of the detector. For the multichannel 
averager, the improvement in signal-to-noise increases as the square root of the number of 
epochs. In either case, the steady-state method provides the amplitude of the response and 
the phase lag of the moment of peak response relative to the onset of the stimulus cycle. 

Measurements of the visually evoked responses at various locations along the midline and 
along a track at right angles to it are illustrated in Figs, la and lb, respectively. The key 
point of interest is that the region of the strongest field was more sharply confined across the 
occipital scalp than were published electric potential patterns, obtained with scalp electrodes. 
This feature is illustrated in Figure lc and d, which show that the region of strong field is 
limited to a span of only 4-5 cm up the midline and extends only 6 cm to the side of it (with 
a suggestion of a second region of strong field appearing farther to the side in Fig. Id). 

It took three more years of work with noisy signals and unreliable sensors for us to realize 
that a neuromagnetic field pattern could be used to determine the three-dimensional location 
of the underlying neuronal activity. This was reported in 1978 in a study of the steady-state 
response of somatosensory cortex to repetitive electrical stimulation of a finger. [2] 

3 Localizing Where the Action Is 

One interesting feature of the visual system is provided by neuronal connections that map 
each solid angle of the visual field onto a specific area of cortex (“retinotopic map”). This 
is only coarsely known in humans. Much more detail, including the existence of multiple 
retinotopic squences, has been established in animals such as monkey and cat through stud- 
ies where electrodes can be placed directly into cortex to observe neuronal activity. The 
retinotopic sequence in humans was discovered in a study of veterans who suffered visual 
defects from wounds. The position of shrapnel in the occipital region of the brain was found 
to have a simple correspondence to the direction in visual space where vision was impaired. [3] 
This relationship was subsequently confirmed during surgical procedures for epileptics when 
electrodes could be placed directly on the pial surface of visual cortex. Electrical responses 
at specific locations across cortex correlated with the placement of a stimulus in the visual 
field. [4] 

Figure 2 illustrates the first non-invasive determination of a retinotopic sequence in 
humans, which were obtained by neuromagnetic techniques. Visual stimuli were presented 
by a grating, in which a periodic array of parallel light, bars alternate with dark ones, and 
the contrast of the pattern is reversed at a fixed rate. By “reversed” we mean that the light 
bars are replaced by dark ones, and vice-versa. To stimulate specific regions of the visual 



Seeing Neuromagnetically 


Page 4 


for measurements of a subject’s reaction time to the appearance of a grating. In other words, 
a behavioral measure shows the same variation with spatial frequency as this physiological 
measure. Taking the value for the neuromagnetic latency and adding about 115 ms predicts 
the reaction time. [7] This suggests a serial processing sequence, in which the occurence of 
peak magnetic field strength is associated with a culminating moment in visual processing. 
Thereafter, a fixed time interval - independent of stimulus properties - elapses before the 
subject releases the button controlling a switch in a reaction time task. We might associate 
the neuromagnetic latency with a sensory processing time and the extra 115 ms as a motor 
response time. 

Why should it take longer for the visual system to respond to a high spatial frequency 
than low, as illustrated by the left-hand part of the data in Fig. 4 ? Breitmeyer [8] was the 
first observed an effect of spatial frequency on reaction time, and he suggested that it may 
well be due to the organization of the visual system as complementary parallel processing 
channels: a ‘sustained’ channel tuned to patterns with high spatial frequency provided that 
movement is not rapid, and a ‘transient’ channel tuned to low spatial frequencies and rapid 
movement, with a more rapid physiological response. This organization had been deduced 
from psychophysical studies of human visual thresholds for gratings of various spatial fre- 
quencies and reversal rates. It has received support from electrophysiological studies on cat 
and monkey as well. The right half of Fig. 4 provides the first physiological evidence for 
humans. [9] 

The hypothesis we tested was the notion that at low contrast reversal rates the ‘sustained’ 
channel dominates the response but at higher rates the ‘transient’ channel would play a more 
important role. This indeed is borne out by the data. [9] For reversal rates above about 20 
Hz, the phase trends display a common slope, indicating that the neuromagnetic latency 
becomes independent of spatial frequency. Moreover, the common slope corresponds to the 
very short latency displayed at lower reversal rates for the lowest spatial frequency. Thus, 
we have evidence that there is a cross-over from pattern- dominated responses at low reversal 
rates to movement-dominated responses at high reversal rates. This is accompanied by an 
increasing processing speed of the visual system. It would be interesting to determine how 
this complementarity in channel performance varies for specific area of visual cortex, by 
using stimuli that are presented in more confined areas of the visual field. 

Okada (unpublished) has shown that the neuromagnetic latency also increases with de- 
crease in luminance or contrast, and reaction times of the same subjects show similar in- 
creases. The results of Okada et al. [10] at low contrast agree with the relation between 
simple reaction time and contrast reported by Harwerth and Levi [11] for evoked potential 
measurements. Therefore, the steady-state latency is a good predictor of reaction time for 
variations of three stimulus characteristics. 

It is interesting that this latency is comparable to the 115 ms we have associated with 
the motor response time. Dividing the total reaction time into more or less equal sensory 
and motor response times makes sense from an operational standpoint. There would be 
comparatively little benefit for nature to invest improved neuronal processing to accelerate 
the speed of a sensory system by, say, a factor of two if there remains a much longer motor 
response time. The return on investment diminishes because the percentage increase in 
overall speed is limited by the slowest system. 
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5 Extravisual Areas of Activity- 

Neurons in medial temporal cortex (MT) of rhesus macaque monkey, which occupies a 
portion of the longitudinal fissure, are known to be selectively sensitive to motion and not to 
orientation of visual patterns. [12] The intriguing notion that it might be possible to detect 
analogous visually- related activity in non- visual areas of the human motivated us to take 
advantage of an extended visit in New York City by Professor Lounasmaa to collaborate on 
a novel study. We sought to develop a stimulus that would permit the detection of neuronal 
activity in response to motion. Since a smoothly moving pattern does not provide a temporal 
reference for detecting an average response that is time locked to a specific stimulus event, we 
sought to detect responses to changes in velocity of a pattern instead. What was found was 
another cortical region of neuronal activity that previously had not been known to respond 
to visual stimulation: an area in the central region, within or near the Rolandic fissure. 

A vertically oriented grating having a spatial frequency of 1.5 cycles per degree of visual 
angle for one subject (3 c/deg for a second subject) was drifted horizontally across an oscil- 
loscope screen, and its speed was sinusoidally modulated at 7 Hz. The modulation produced 
a “velocity contrast” of 50%, expressed by the modulation amplitude as a percentage of the 
mean velocity. With the subject fixated on the center of the screen, the modulation was 
keep sufficiently low that the grating drifted uni directionally. Neuromagnetic responses were 
observed over a wide area of the scalp by Lounasmaa et al., [13] as illustrated in the sagittal 
view of Fig. 5a. Two areas of strong field, each having a common phase lag but differing in 
phase by 180 deg, appear over the lower occipital region. Only the right hemisphere pattern 
appears in the figure. A weaker region of apparently random phase appears just above. In 
separate studies of half-field presentations this lower and upper area can be identified as the 
regions of inward and outward field from activity in visual cortex of the right hemisphere. A 
similar pattern appears on the left hemisphere, but with field directions reversed. If each is 
modeled by a current dipole source, the observed pattern is explained by having the dipoles 
tipped slightly from the horizontal so that the two strong field extrema over the upper oc- 
cipital area largely superimpose. This would account for the weak magnetic field in that 
region for full visual-field present aton. 

When the field patterns of the current dipoles representing left and right occipital sources 
are subtracted from the measured pattern (Fig. 5b), an interesting feature remains in the 
data recorded over the Rolandic fissure: there are two regions within which the phase is 
consistent but with a 180 deg difference between the regions. One region lies near the border 
of the diagram just posterior to the vertex (phasors pointing toward the lower left) and the 
other lies to the lower right of it (phasors pointing toward the upper right). A similar pattern 
is found over the left hemisphere (not shown), except the field directions relative to the scalp 
are reversed. Thus each source in the left and right central areas may be modeled by an 
equivalent current dipole lying under the midpoint of the pattern. The corresponding phaser 
pattern being shown in Fig. 5c. When this predicted field pattern is subtracted Rom the data 
of Fig. 5b, the resulting residual field pattern in Fig. 5d is without any regions of appreciable 
amplitude and consistent phase, implying that no additional sources of appreciable strength 
contribute to it. 

The central sources in the right hemisphere lies within or close to the Rolandic fissure, 
whose locus is well established for the subject characterized in Fig. 5. The source lies 2 cm 
lower along the Rolandic fissure than the projection for motor and somatosensory response 
for voluntary ballistic flexure of the left index finger. [14] Its orientation is perpendicular 
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to the fissure, and its depth is 3 cm beneath the scalp, thus placing it within the fissure. 
To our knowledge, this was the first identification of visually related activity in this area. 
An interesting feature of these central magnetic fields is that activity appears whether the 
stimulus is presented in right or left halves of the visual field. Such bilateral responses 
distinguish this neuronal activity in the central region from activity in visual cortex. 

Thus neuromagnetic measurements made it possible to detect a new region of activity 
in an extra- visual area that had not been identified previously by microelectrode studies on 
primates. This activity is evoked by changes in the velocity of the stimulus, which would seem 
to imply that the neuronal population is involved in a tracking function. There is preliminary 
evidence that the latency of this activity is comparatively long (approximately 200 ms), and 
we have speculated elsewhere that its function may well be to follow the processing of images 
in visual cortex for comparison purposes, [13] such as determining whether an image moving 
across the retina is due to movement of the object or of the eye. Thus, proximity of this 
visually active region to motor cortex suggests that it may be the site to which reafferent 
signals are distributed for such comparison purposes. 

6 Visual Imagery 

One of our principal interests is to extend the capability of magnetic studies to higher levels 
of brain function, which includes tasks such as memory, cognition, and making decisions. 
A first step was the study of the so-called “P-300” signal, which is commonly thought to 
indicate an updating of a person’s expectations for the immediate future.[15] A source for 
this, as well as for a related “N-200” signal was found to be the hippocampal formation, deep 
within temporal lobe. However, our emphasis in this article is on the visual system, so we 
shall direct our attention to a new topic we have recently begun to explore: mental imagery. 
Simply put, the question is whether we use much the same neural machinery in mentally 
visualizing the face of someone familiar as we do in real-time vision. Is visual cortex involved 
in mental imagery? 

We expected that cortical activity in such a task would be no stronger than for normal 
visual processing, so the experimental difficulty would be in detecting such a weak magnetic 
signal when there was no way to exactly repeat it and use synchronous averaging. Instead, 
a “lever” was used whereby the much stronger alpha rhythm (8-13 Hz) was monitored 
magnetically and changes in its spatial distribution during specific tasks were noted. This 
approach was motivated in part by previous studies of the scalp potential, [16] [17] where 
the power in different bandwidths is differentially attenuated across the scalp depending on 
specific cognitive tasks that a subject is performing. 

The task we chose is one where the subject sees a sequence of three or more abstract 
figures, each shown for 1 sec, and separated by 0.3 sec. Following that by a few seconds 
a “probe” figure is shown for only 0.1 sec. The instructions are to press one reaction time 
button if the probe was seen previously or a second button if not. The magnetic field is 
recorded through the entire sequence. Then the waveform representing the variance about 
the average magnetic signal within the bandwidth 8-12 Hz is computed for each epoch. This 
variance represents power in the the alpha bandwidth that is not time-locked to the visual 
presentations. The variance waveforms were averaged across 30 epochs and smoothed by a 
low-pass filter at 8 Hz. , 

Figure 6 illustrates two examples of the time-course of this alpha power.[18] One (broken 
line) shows how alpha power over the occipital area of the scalp is modulated when a subject 
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field the grating was shown only in the right visual field, and it was contoured as the right 
half of an annulus. This produced a stimulus with a uniform angular displacement from the 
center of the visual field, when the subject fixated his gaze at a small point of light at the 
center of the corresponding annulus. 

The lower portion of Fig. 2 shows the steady-state responses as phasors, where the center 
of each vector indicates a measurement position, its length gives the amplitude of the re- 
sponse, and its orientation (measured clockwise from the right horizontal direction) indicates 
the phase lag. For each, large ovals are drawn around the two areas of strongest field. The 
field normal to the scalp reverses direction from one to the other (phasors are antiparallel). 
The nearly vertical alignment of these “field extrema” indicates the corresponding neuronal 
source lies between them, with current oriented nearly horizontal. This is expected on the 
basis of the cruciform model of cortex, if the local neuronal current producing the field flows 
perpendicular to the cortical surface. Then the field of the cortical layer forming the floor 
of the calcarine fissure in the left hemisphere is largely cancelled by the field of the cortical 
layer forming the ceiling just above. Only cortex forming the left wall of the longitudinal 
fissure provides a net contribution, and its current would be horizontal. The neuronal source 
of the observed field is believed to be the pyramidal cells within cortex. They are the only 
population that has a preferred orientation for their dendrites, and that alignment is indeed 
perpendicular to the cortical surface. 

If the neuronal source is modeled as a small area of activity (a “current dipole”), the 
distance between the extrema in comparison with the radius of the head predicts the depth 
of the dipole beneath the scalp. [5] Of particular interest is the fact that the span between 
extrema increases monotonically as the stimulus is presented more peripherally. This im- 
plies that the corresponding neuronal source lies deeper within the brain, as illustrated in 
Fig. 3. [6] While these first measurements would be considered comparatively crude by to- 
day’s standards for magnetic localization, the sequence of positions is in general agreement 
with the expected locus of visual cortex in the sagittal plane, and with more recent studies 
by invasive techniques. 

4 Functional Properties: Pattern versus Movement 

Steady-state investigations of sensory processes have not been as popular as transient mea- 
surements. However, the ability to manipulate the rate of stimulation over a large range 
makes the steady-state approach attractive for exploring temporal properties of neuronal 
mechanisms. One example is illustrated by the data in Fig. 4. Here the phase lag for peak 
field relative to the stimulus is plotted versus contrast reversal rate. Results are shown 
for a large contrast-reversing grating subtending an angle of 10 deg horizontally and 5 deg 
vertically. For low reversal rates, phase lag increases linearly with reversal rate. Moreover, 
the slope increases with the spatial frequency of the grating. The linear relationship can be 
ascribed to a fixed delay between stimulus and peak response. When stimuli are presented 
with shorter intervals, the response following the first stimulus lies farther back in the stim- 
ulus cycle, hence a greater phase lag. Values of phase lag greater than 2 tt indicate that 
a second stimulus has been applied before the peak response to the first one takes place. 
Values greater than 47r indicate that two stimuli are interposed before the response to a 
previous one. * 

Is the increase of response latency with increasing spatial frequency shown in Fig. 4 at 
low reversal rates meaningful? The answer is affirmative, for exactly this variation is found 
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is instructed to press a button when the probe appears, whether or not the figure was 
previously seen. Alpha power is comparatively high during the second preceding the probe’s 
appearance (at the time origin) but there is a brief suppression lasting about half a second 
following the appearance. The reaction time for the button press in this ‘simple reaction 
time’ task is indicated by an arrow, with the standard deviation given by a horizontal bar. 
In this case the reaction time is about 0.3 sec. Quite a different behavior is seen when the 
subject mentally compares the probe with figures previously seen (solid line). Here again 
there is a prompt suppression of alpha power when the probe appears, but suppression is 
maintained for more than 1 sec, lasting until about the time a reaction time button is pressed 
to indicate a match or no match in this ‘choice reaction time’ task. Thus, the duration of 
suppression during this task of mental imagery is comparable to the time subjects take to 
make a decision. 

What is particularly interesting is the spatial distribution of the alpha suppression across 
the posterior area of the scalp. Figure 7 illustrates this for three subjects, where positions 
across the scalp are projected onto a plane for representing the power patterns. Fig. 7b shows 
the distribution of “baseline alpha” power prior to presentation of the probe. Generally 
there are two regions of strongest field, over left and right hemispheres. This is not so well 
represented by subject BS, and not at all by subject CS whose alpha power is two orders of 
magnitude lower. These differences across subjects illustrate in interesting but unexplained 
property of the alpha rhythm: Its strength varies greatly across individuals. During the 
mental imagery task, there is general suppression over the entire area, although some spatial 
structure in the power distribution is still evident in Fig. 7c. This is called the residual 
alpha” distribution. 

What is particular relevant in our study is the relative suppression of alpha, expressed 
as the distribution of “relative residual alpha” in Fig. 7d. This is defined as the ratio of 
the residual to baseline alpha power. If there were uniform suppression across the scalp, 
the relative suppression would show no spatial features. On the contrary, Fig. 7d shows for 
subjects LK and BS that there is a strong spatial dependence of the relative suppression, 
with the most pronounced effect being near the midline of the occipital scalp. (Subject CS 
has such weak baseline alpha power that there is very little suppression anywhere, and so we 
do not consider his case further.) The conclusion we draw is that the neuronal sources that 
would produce a magnetic field in this midline region of enhanced suppression lie within the 
area of visual cortex. Therefore, the visual cortex is indeed involved in the process of mental 
imagery. 

It remains to be determined whether this activity is related to retrieval of image infor- 
mation from short-term memory associated with this cortical area, or whether perhaps it 
is related to the subject’s forming his perception. In either case, the use of alpha rhythm 
as a leverage device opens many prospects for more specifically targeted studies of mental 
imagery. The technique is also applicable to studies of cognitive tasks in other modalities. 


7 Retrospective 

While this was not intended as a complete survey of neuromagnetic studies of visual pro- 
cesses, it is quite representative. In fact, given the precision of results obtained using these 
methods in other sensory modalities, it is somewhat surprising that more has not been ac- 
complished over the past dozen years. We would attribute this to several factors. Previously 
it has been difficult to provide visual stimulaiton without producing magnetic artifacts. Also 
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it was difficult to keep a subject in a fixed position while measuring the field over a broad 
region of the occiptial, parietal, and temporal areas. The canted sensors in the 7- sensor and 
24-sensor systems developed at the Low Temperature Laboratory of the Helsinki University 
of Technology considerably ameliorate this problem. Other recent improvements include 
CryoSQUID sensors [19] that can be oriented in any direction, video projectors that can be 
used outside a magnetically shielded room to project the image inside, and versatile gantries 
that ease the task of positioning the dewar containing the magnetic sensors. These advances 
in instrumentation promise to markedly enhance the productivity of neuromagnetic studies 
of the human visual system and encourage more researchers to become involved. 

In this brief presentation we have concentrated on aspects of visual studies that we have 
found exciting and which hold promise for many more fruitful investigations by magnetic 
techniques. The combination of a physical technique nutured in the laboratories of low 
temperature physicists and psychological techniques for probing brain functions is powerful. 
It is with our warm regards to a colleague, friend, and fine competetor that we dedicate this 
article to Professor Olli V. Lounasmaa. 
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9 Figure Captions 

Fig. 1. Time variation of the magnetic field evoked by a bar grating flashed on for 50 ms 
and off for 50 ms to produce a 10-Hz flicker. Magnetic response (A) and response power (C) 
at various distances along the midline above the inion and the response (B) and response 
power (D) as a function of distance to the left of the midline. The archaic unit of the gauss, 
shown as a unit of measure beside the calibration bar for field strength, is equivalent to 10 -4 
tesla. (Brenner et al. 1975, with permission of American Association for the Advancement 
of Science) 

Fig. 2. Visual stimuli presented at increasing distance from the subject’s fixation point 
(dot) and the corresponding evoked field in the regions of strongest responses, represented 
by phasors. The length of a vector indicates the amplitude of response at that location, and 
its angle measured counterclockwise from the right horizontal axis indicates the phase lag. 
(adapted from Maclin et al. 1983) 

Fig. 3. SaggitaL view of a representative head (anterior to the left), displaying the three 
positions deduced for the neuronal sources that account for the field patterns of Fig. 2. Star 
indicates the center of curvature for the spherical head used in fitting a current dipole source 
model to the data. Dashed line indicates the boundary of visual cortex according to a stan- 
dard anatomical atlas. 

Fig. 4. Trends for the phase lag of the steady-state neuronal response to contrast reversal 
stimuli at various repetition rates, for a grating of the indicated spatial frequency. Phases 
for rates above 20 Hz are decreased by 2 tt in order to better show the common trend, (from 
Williamson et al. 1979) 

Fig. 5. (a) Sagittal views of a spherical representation of the subject’s scalp, with phasors 
characterizing the response to a velocity-modulated moving grating. Strong activity is ev- 
ident in the occipital area (lower left) and weaker activity with consistent phase lags are 
in the central area (upper right), (b) Resulting pattern when phasors corresponding to a 
current dipole in each of the left and right occipital areas are subtracted, (c) Field pattern 
for a current dipole placed in the Rolandic fissure, showing the inward and outward field 
regions by phasors having reversed directions, (d) Residual field pattern when the field in 
(c) is subtracted from the data in (b). No preferred phase is evident for the phasors in the 
field region of the Rolandic sulcus. Symbols: square indicates inion; circle, right ear canal; 
triangle, vertex. (Lounasmaa et al. 1985) 

Fig. 6. Amplitude of alpha power near the midline of the occipital scalp is suppressed mo- 
mentarily when an abstract figure is shown (broken line) but suppression is prolonged when 
the subject mentally compares the figure' with three previously viewed (solid line). Reaction 
times to indicate task completion also differ: a quick response for the former (Simple RT) 
but much longer for the latter (Choice RT). (Kaufman et al. 1990) 

Fig. 7. (a) Azimuthal equal distance projection onto a plane, in which distances in any di- 
rection from a point on a spherical scalp are represented by equal distances in corresponding 
directions from the overlying point on a tangential plane, (b) Distribution of average base- 
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line alpha for such a projection, with z-axis of this plane corresponding to a horizontal line 
through the inion while the y axis is the projection of the midline. Baseline alpha is defined 
as the average alpha power observed within a 200 msec interval 100 msec prior to presenta- 
tion of the visual probe, (c) Distribution of residual alpha , defined as power averaged over a 
100 msec interval centered on the moment of maximum suppression, (d) Distribution of the 
ratio of residual to baseline alpha, which defines the relative residual alpha. (Adapted from 
Kaufman et al. 1990) 



